You are here: R workshop >> Session 6 Mapping >> mapping demo


Nordic Oikos 2018 - R workshop - Session 6

Session 6 focuses on working with environment layers, mapping, cropping and masking layers using the Raster R-package and other tools.


GBIF data for taxon liverleaf (blåveis:no) - from Norway

library('rgbif') # rOpenSci r-package for GBIF data
library('mapr') # rOpenSci r-package for mapping (occurrence data)
sp_name <- "Hepatica nobilis"; kingdom <- "Plantae" # liverleaf (blaaveis:no), taxonKey=5371699
key <- name_backbone(name=sp_name, kingdom=kingdom)$speciesKey
sp <- occ_search(taxonKey=key, return="data", hasCoordinate=TRUE, country="NO", limit=100)
sp_m <- sp[c("name", "catalogNumber", "decimalLongitude","decimalLatitude", "basisOfRecord", "year", "municipality", "taxonKey", "occurrenceID")] ## Subset columns (for more useful map pop-up)
#gbifmap(sp, region = "norway")
map_leaflet(sp_m, "decimalLongitude", "decimalLatitude", size=2, color="blue")
Map GBIF data for Hepatica nobilis from Norway

Map GBIF data for Hepatica nobilis from Norway


Extract coordinates suitable for e.g. Maxent

xy <- sp[c("decimalLongitude","decimalLatitude")] ## Extract only the coordinates
sp_xy <- sp[c("species", "decimalLongitude","decimalLatitude")] ## Input format for Maxent
#head(sp_xy, n=5) ## preview first 5 records

Write dataframe to file (useful for Maxent etc.)

write.table(sp_xy, file="./demo_data/sp_xy.txt", sep="\t", row.names=FALSE, qmethod="double") ## for Maxent
#readLines("./demo_data/sp_xy.txt", n=10)
write.table(sp, file="./demo_data/sp.txt", sep="\t", row.names=FALSE, qmethod="double") ## dataframe
#readLines("./demo_data/sp.txt", n=10)

Read data file back into R

sp_xy <- read.delim("./demo_data/sp_xy.txt", header=TRUE, dec=".", stringsAsFactors=FALSE)
sp <- read.delim("./demo_data/sp.txt", header=TRUE, dec=".", stringsAsFactors=FALSE) ## dataframe

Get (detailed) administrative borders for Norway from GADM

Very slow to plot, many details…

library(raster)
gadm_norway_1 <- getData('GADM', country='NOR', level=1, path="./demo_data") ## level 0,1,2,...
plot(gadm_norway_1, main="Adm. Boundaries Norway Level 1")
points(xy, col='blue', pch=20) ## plot species occurrence points to the map (smaller dots)
#legend("bottomright", title = "Legend", legend = "H. nobilis", pch = 20, pt.bg = "blue", bty = "n")
GADM Norway

GADM Norway


Get (simpler) country borders for Norway from maptools

library(maptools)
library(rgeos)
data(wrld_simpl) ## vector/factor of ISO2, ISO2, NAME, REGION, SUBREGION, LON, LAT, ...
norway_mask <- subset(wrld_simpl, NAME=="Norway") ## extract the (simpler) border for Norway
plot(norway_mask, axes=FALSE, border="#777777") ## 
points(xy, col='blue', pch=20, cex=1) # plot species occurrence points to the map
title("Country mask for Norway from wrld_simpl")
legend("bottomright", title = "Legend", legend = "Occurrences", pch = 20, col="blue", cex = 0.9)
Border for Norway from maptool:wrld_simpl

Border for Norway from maptool:wrld_simpl


Admin borders from DIVA-GIS

DIVA-GIS data by country: [http://www.diva-gis.org/gdata]. DIVA-GIS admin borders for Norway: [http://biogeo.ucdavis.edu/data/diva/adm/NOR_adm.zip]. Notice that these are the same (detailed) borders as available from GADM.

## Download DIVA-GIS admin borders for Norway
if (!file.exists('./demo_data/NOR_adm.zip')) {
    download.file('http://biogeo.ucdavis.edu/data/diva/adm/NOR_adm.zip', './demo_data/NOR_adm.zip')
    dir.create(file.path("./demo_data/NOR_adm")) ## create a folder for shape files
    unzip('./demo_data/NOR_adm.zip', exdir='./demo_data/NOR_adm')
}
#NOR_adm1 <- shapefile('./demo_data/NOR_adm/NOR_adm1.shp') ## county borders
#plot(NOR_adm1, axes=FALSE, border="#666666")
DIVA-GIS NOR_adm1 (source: GADM)

DIVA-GIS NOR_adm1 (source: GADM)

library(raster)
library(sp)
#NOR_adm0 <- shapefile('./demo_data/NOR_adm/NOR_adm0.shp') ## country border
#NOR_adm1 <- shapefile('./demo_data/NOR_adm/NOR_adm1.shp') ## county borders
#NOR_adm2 <- shapefile('./demo_data/NOR_adm/NOR_adm2.shp') ## municipality borders
par(mfrow=c(1,3)) ## combining plots: nrows, ncols
##
plot(NOR_adm0, axes=FALSE, border="#666666")
points(xy, col='blue', pch=20, cex=1) ## pecies occurrence
title("NOR_adm0, country")
##
plot(NOR_adm1, axes=FALSE, border="#666666") ## may set line width -- lwd=INT, but increases rendering time substatially
points(xy, col='blue', pch=20, cex=1) ## species occurrence
title("NOR_adm1, county")
##
plot(NOR_adm2, axes=FALSE, border="#666666")
points(xy, col='blue', pch=20, cex=1) ## species occurrence
title("NOR_adm2, municipality")
#legend("bottomright", title = "Legend", legend = "Occurrences", pch = 20, col="blue", cex = 0.9)
##
DIVA-GIS NOR_adm (source: GADM)

DIVA-GIS NOR_adm (source: GADM)


Read environment layer from WorldClim into R

For more information about WorldClim see session 5.

require(raster) # spatial raster data management, works well with dismo
env <- getData('worldclim', var='bio', res=10) # 10 degree grid (approx 18.5 km, 342 km2 at equator) 85 MByte

Plot environment layers and species occurrences on a map

plot(env, 1, main=NULL, axes=FALSE) ## could add title here with main="Title"
title(main = bquote(italic(.(sp_name)) ~occurrences~on~Annual~mean~temperature~'(dCx10)'))
#plot(gadm_norway, add=TRUE) ## add admin county borders
points(xy, col='blue', pch=20) # plot species occurrence points to the map (smaller dots)
Bioclim 1, Annual mean temperature

Bioclim 1, Annual mean temperature

plot(env, 12, main=NULL, axes=FALSE) # plot bioclim layer, BIO12 = Annual Precipitation
title(main = bquote(italic(.(sp_name)) ~occurrences~plotted~on~Annual~precipitation~'(mm)'))
axis(side=2, tck = -0.04, col.ticks="gray") ## add axis only left
points(xy, col='blue') # plot species occurrence points to the map
Bioclim 12, Annual precepitation

Bioclim 12, Annual precepitation



Crop and mask raster layers

Cutting large (global) environment layers to a smaller extent can save significant memory. If your species occurrence data are limited to a region (e.g. Norway, Scandinavia or similar) you might reduce computation time significantly by cropping your environment layers appropriatly.

Cut environment layer(s) to extent (result is always a square)

library(raster)
library(rgdal)
ext <- extent(3,35,54,72) ## minLon=3, maxLon=35, minLat=54, MaxLat=72 for Scandinavia
env_cut <- crop(env, ext, snap="out") ## square output
#env_cut <- crop(env, ext, snap="out", filename="./demo_data/env_cut.tif") ## add filename to save result
plot(env_cut)
Bioclim layers cropped to Scandinavia

Bioclim layers cropped to Scandinavia

Cut environment layer(s) to a mask (from a shapefile or other vector data)

data(wrld_simpl) ## here you can also read in your OWN vector data (e.g. study area)
norway_mask <- subset(wrld_simpl, NAME=="Norway")
env_crop <- crop(env, extent(norway_mask)) ## crop gives a square (cut to the extent of the mask)
env_mask <- mask(env_crop, norway_mask) ## mask removes data outside the mask
plot(env_mask)
plot(norway_mask, add=TRUE, lwd=2)
Bioclim layers masked to Norway

Bioclim layers masked to Norway


Plot with extent Scandinavia (using zoom)

Using zoom, the raster data in R workspace environment is still the same size. You only zoom into the region of interest for more useful maps.

Bioclim 12, Annual precepitation

Bioclim 12, Annual precepitation


Map of the cropped and masked raster layers

par(mfrow=c(1,2)) ## combining two plot with par(n_rows, n_columns)
##plot(env_cut$bio12); title(main="Bio12 cut to Scandinavia")
plot(env_crop$bio12); title(main="Bio12 cropped to Norway")
points(xy, col='blue', pch=20) ## add species occurrence points to cropped map
plot(env_mask$bio12); title(main="Bio12 masked to Norway")
points(xy, col='blue', pch=20) ## add species occurrence points to masked map
Bioclim 12 (Annual precepitation) respectively cropped and masked to Norway

Bioclim 12 (Annual precepitation) respectively cropped and masked to Norway



Testing remote sensing image for Trondheim downloaded from Landsat

I used the USGS LandsatLook Viewer and Sentinel2Look Viewer to download sattelite data for Trondheim.

rs_l <- brick('./demo_data/landsat_trondheim_web_mercartor_wgs84.tif')
rs_s <- brick('./demo_data/sentinel_trondheim_web_mercartor_wgs84.tif')
nlayers(rs_l); nlayers(rs_s) ## 3 layers
crs(rs_l); crs(rs_s) ## +proj=merc +a=6378137 +b=6378137 ... +units=m ...
ncell(rs_l); ncell(rs_s) ## rs_l = 100738 ## rs_s = 204530
dim(rs_l); dim(rs_s) ## rs_l = 241 418   3 ## rs_s = 362 565   3
res(rs_l);res(rs_s) ## rs_l = 30.09113 30.09113 ## rs_s = 20.03305 20.03305
#par(mfrow=c(2,1)) ## combining two plot with par(n_rows, n_columns)
#plotRGB(rs_l, stretch="lin", axes=FALSE, main="Landsat True Color Composite")
plotRGB(rs_s, stretch="lin", axes=FALSE, main="Sentinel True Color Composite")
Remote sensing data, sentinel, Trondheim

Remote sensing data, sentinel, Trondheim


Base map from Google

Base map from Google, with H. nobilis, Trondheim, Oslo

Base map from Google, with H. nobilis, Trondheim, Oslo


Size of environment layer can be LARGE if using the finer resolutions

## object.size(env) ## read the space allocated in memory for an environment variable
## format(object.size(env), units = "auto") ## Auto reports multiples of 1024
## format(object.size(env), units = "auto", standard = "SI") ## SI use multiples of 1000
cat("Size of env =", format(object.size(env), units = "auto"))
cat("\nSize of env_cut =", format(object.size(env_cut), units = "auto"))
cat("\nSize of env_mask =", format(object.size(env_mask), units = "auto"))
cat("\nSize of gadm_norway_1 =", format(object.size(gadm_norway_1), units = "auto"))
#rm(env) ## save memory - especially useful if using finer resolutions
  • Size of env = 235.5 Kb
  • Size of env_cut = 13.7 Kb
  • Size of env_mask = 940.8 Kb
  • Size of gadm_norway_1 = 12.8 Mb


The BioClim layers:

  • BIO1 = Annual Mean Temperature
  • BIO2 = Mean Diurnal Range (Mean of monthly (max temp – min temp))
  • BIO3 = Isothermality (BIO2/BIO7) (* 100)
  • BIO4 = Temperature Seasonality (standard deviation *100)
  • BIO5 = Max Temperature of Warmest Month
  • BIO6 = Min Temperature of Coldest Month
  • BIO7 = Temperature Annual Range (BIO5-BIO6)
  • BIO8 = Mean Temperature of Wettest Quarter
  • BIO9 = Mean Temperature of Driest Quarter
  • BIO10 = Mean Temperature of Warmest Quarter
  • BIO11 = Mean Temperature of Coldest Quarter
  • BIO12 = Annual Precipitation
  • BIO13 = Precipitation of Wettest Month
  • BIO14 = Precipitation of Driest Month
  • BIO15 = Precipitation Seasonality (Coe cient of Variation)
  • BIO16 = Precipitation of Wettest Quarter
  • BIO17 = Precipitation of Driest Quarter
  • BIO18 = Precipitation of Warmest Quarter
  • BIO19 = Precipitation of Coldest Quarter

GBIF data for taxon liverleaf (blåveis:no) - from Trondheim

library('rgbif') # rOpenSci r-package for GBIF data
library('mapr') # rOpenSci r-package for mapping (occurrence data)
sp_name <- "Hepatica nobilis"; kingdom <- "Plantae" # liverleaf (blaaveis:no), taxonKey=5371699
key <- name_backbone(name=sp_name, kingdom=kingdom)$speciesKey
bb <- c(10.2,63.3,10.6,63.5) # Trondheim
#bb <- c(5.25, 60.3, 5.4, 60.4) # Bergen
#bb <- c(18.7, 69.6, 19.2, 69.8) # Tromsoe
#bb <- c(10.6, 59.9, 10.9, 60.0) # Oslo
sp_bb <- occ_search(taxonKey=key, return="data", hasCoordinate=TRUE, country="NO", geometry=bb, limit=100)
sp_bb_m <- sp_bb[c("name", "catalogNumber", "decimalLongitude","decimalLatitude", "basisOfRecord", "year", "municipality", "taxonKey", "occurrenceID")] ## Subset columns
map_leaflet(sp_bb_m, "decimalLongitude", "decimalLatitude", size=2, color="blue")
Map GBIF data with bounding box for Trondheim

Map GBIF data with bounding box for Trondheim


Navigate back to GitHub project home or GitHub.io html pages.